Linear Regressions with Fraction Diversity

# Free-Living samples only 
div_measures_free <- otu_alphadiv %>%
  dplyr::select(fraction, mean, measure, frac_bacprod, fracprod_per_cell_noinf) %>%
  filter(fraction == "Free") %>%
  dplyr::select(-fraction)
# Particle-associated samples only 
div_measures_part <- otu_alphadiv %>%
  dplyr::select(fraction, mean, measure, frac_bacprod, fracprod_per_cell_noinf) %>%
  filter(fraction == "Particle") %>%
  dplyr::select(-fraction)

## All of the samples!
div_measures_all <- otu_alphadiv %>%
  dplyr::select(mean, measure, frac_bacprod, fracprod_per_cell_noinf)

########################################################
##############  Particle: Community-wide 
lm_divs_particle_comm <- div_measures_part %>%
  dplyr::select(-fracprod_per_cell_noinf) %>%
  group_by(measure) %>% 
  do(glance(lm(frac_bacprod ~ mean, data = .))) %>% 
  mutate(fraction = "Particle", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  Particle: Community-wide 
lm_divs_particle_percap <- div_measures_part %>%
  dplyr::select(-frac_bacprod) %>%
  group_by(measure) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ mean, data = .))) %>% 
  mutate(fraction = "Particle", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

########################################################
##############  Free: Community-wide 
lm_divs_free_comm <- div_measures_free %>%
  dplyr::select(-fracprod_per_cell_noinf) %>%
  group_by(measure) %>% 
  do(glance(lm(frac_bacprod ~ mean, data = .))) %>% 
  mutate(fraction = "Free", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  Free: Community-wide 
lm_divs_free_percap <- div_measures_free %>%
  dplyr::select(-frac_bacprod) %>%
  group_by(measure) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ mean, data = .))) %>% 
  mutate(fraction = "Free", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

########################################################
##############  All Samples: Community-wide 
lm_divs_all_comm <- div_measures_all %>%
  dplyr::select(-fracprod_per_cell_noinf) %>%
  group_by(measure) %>% 
  do(glance(lm(frac_bacprod ~ mean, data = .))) %>% 
  mutate(fraction = "All Samples", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  All Samples: Community-wide 
lm_divs_all_percap <- div_measures_all %>%
  dplyr::select(-frac_bacprod) %>%
  group_by(measure) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ mean, data = .))) %>% 
  mutate(fraction = "All Samples", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, measure, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 


## Filtered PCA Linear Model Results
# Put all of the PCA and environmental  linear model results together into one dataframe 
lm_div_results <- 
  bind_rows(lm_divs_particle_comm, lm_divs_particle_percap,
          lm_divs_free_comm, lm_divs_free_percap,     
          lm_divs_all_comm, lm_divs_all_percap) %>% 
  dplyr::rename(independent_var = measure,
                dependent_var = dependent) 

See supplemental table 1 header for the table output!

Cross Validation

Cross validation will provide a better estimate for the adjusted R-squared value of our diversity regressions.

#################################### Community-Wide Production #################################### 
################### Richness ################### 
lm_prod_vs_rich_PA <- lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Particle"))

summary(lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(richness, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1271 -2.8865 -0.4649  3.7537  9.8401 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -11.269531   5.717450  -1.971  0.07701 . 
## mean          0.038125   0.009868   3.863  0.00314 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.476 on 10 degrees of freedom
## Multiple R-squared:  0.5988, Adjusted R-squared:  0.5587 
## F-statistic: 14.93 on 1 and 10 DF,  p-value: 0.003143

summary(lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(richness, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.094 -11.112  -2.755   8.611  33.308 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.40793   21.62006   0.250    0.808
## mean         0.05511    0.06207   0.888    0.395
## 
## Residual standard error: 17.7 on 10 degrees of freedom
## Multiple R-squared:  0.07306,    Adjusted R-squared:  -0.01963 
## F-statistic: 0.7882 on 1 and 10 DF,  p-value: 0.3955

# Without the 2 high samples 
summary(lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Particle" & mean < 700)))

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(richness, fraction == 
##     "Particle" & mean < 700))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7578 -1.6875 -0.6392  2.8807  6.3336 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.44152    8.50828  -0.640    0.540
## mean         0.02559    0.01701   1.505    0.171
## 
## Residual standard error: 4.36 on 8 degrees of freedom
## Multiple R-squared:  0.2205, Adjusted R-squared:  0.1231 
## F-statistic: 2.264 on 1 and 8 DF,  p-value: 0.1709

rich_test_df <- 
  richness %>%
  mutate(frac_boolean = (fraction == "Particle")*1) 

# Correlation of richness between PA and FL
rich_test_df_sorted <- arrange(rich_test_df, norep_water_name)
rich_test_df_sorted_PA <- filter(rich_test_df_sorted, fraction == "Particle")
rich_test_df_sorted_FL <- filter(rich_test_df_sorted, fraction == "Free")
cor.test(rich_test_df_sorted_PA$mean, rich_test_df_sorted_FL$mean)

## 
##  Pearson's product-moment correlation
## 
## data:  rich_test_df_sorted_PA$mean and rich_test_df_sorted_FL$mean
## t = 0.79099, df = 10, p-value = 0.4473
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3848365  0.7167446
## sample estimates:
##       cor 
## 0.2426582

# Look into car::vif() for multicollinearity 
cor(rich_test_df$mean, rich_test_df$frac_boolean)

## [1] 0.6508802

summary(lm(frac_bacprod ~ mean + frac_boolean, data = rich_test_df))

## 
## Call:
## lm(formula = frac_bacprod ~ mean + frac_boolean, data = rich_test_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.525  -5.808  -1.127   5.678  34.828 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    9.95501    7.87833   1.264  0.22022   
## mean           0.04167    0.02055   2.028  0.05545 . 
## frac_boolean -23.20045    6.89481  -3.365  0.00293 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.82 on 21 degrees of freedom
## Multiple R-squared:  0.3512, Adjusted R-squared:  0.2894 
## F-statistic: 5.684 on 2 and 21 DF,  p-value: 0.01064

PA_rich_residplot <- plot_residuals(lm_prod_vs_rich_PA, filter(richness, fraction == "Particle")$frac_bacprod, "Community Production vs Particle-Associated Richness")

## [1] "There are no high leverage points in this model."

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_prod_vs_rich_PA <- train(
      frac_bacprod ~ mean, data = filter(richness, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_prod_vs_rich_PA$results   # Particle Community-Wide CV results 

##   intercept     RMSE  Rsquared    MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.878977 0.6666884 5.5294 2.263403  0.2802459 1.908415

################### Shannon Entropy ################### 
lm_prod_vs_shannon_PA <- lm(frac_bacprod ~ mean, data = filter(shannon, fraction == "Particle"))
cv_lm_prod_vs_shannon_PA <- train(
      frac_bacprod ~ mean, data = filter(shannon, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_prod_vs_shannon_PA$results   # Particle Community-Wide CV results 

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.620109 0.6805392 5.226135 2.393256   0.285612 2.069473

################### Inverse Simpson ################### 
lm_prod_vs_invsimps_PA <- lm(frac_bacprod ~ mean, data = filter(invsimps, fraction == "Particle"))
summary(lm_prod_vs_invsimps_PA)

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(invsimps, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5589 -2.1093 -0.1969  0.8752  7.8282 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.45199    2.43866  -0.185  0.85666    
## mean         0.29344    0.05781   5.076  0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.571 on 10 degrees of freedom
## Multiple R-squared:  0.7204, Adjusted R-squared:  0.6925 
## F-statistic: 25.77 on 1 and 10 DF,  p-value: 0.0004804

# Without the 2 high samples 
summary(lm(frac_bacprod ~ mean, data = filter(invsimps, fraction == "Particle" & mean < 70)))

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(invsimps, fraction == 
##     "Particle" & mean < 70))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3309 -1.6323 -0.3605  1.4988  6.5349 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  1.61837    2.49960   0.647   0.5355  
## mean         0.20554    0.08156   2.520   0.0358 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.688 on 8 degrees of freedom
## Multiple R-squared:  0.4425, Adjusted R-squared:  0.3729 
## F-statistic: 6.351 on 1 and 8 DF,  p-value: 0.03581

PA_invsimps_residplot <- plot_residuals(lm_prod_vs_invsimps_PA, filter(invsimps, fraction == "Particle")$frac_bacprod, "Community Production vs Particle-Associated Inverse Simpson")

## [1] "There are no high leverage points in this model."

invsimps_test_df <- 
  invsimps %>%
  mutate(frac_boolean = (fraction == "Particle")*1) 

cor(invsimps_test_df$mean, invsimps_test_df$frac_boolean)

## [1] 0.3165699

summary(lm(frac_bacprod ~ mean + frac_boolean, data = invsimps_test_df))

## 
## Call:
## lm(formula = frac_bacprod ~ mean + frac_boolean, data = invsimps_test_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.355  -5.937  -0.176   2.676  37.452 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   16.3742     5.2053   3.146  0.00488 **
## mean           0.3189     0.1527   2.089  0.04907 * 
## frac_boolean -17.7310     5.4908  -3.229  0.00402 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.76 on 21 degrees of freedom
## Multiple R-squared:  0.3576, Adjusted R-squared:  0.2965 
## F-statistic: 5.846 on 2 and 21 DF,  p-value: 0.009586

# Correlation of inverse Simpson's index between PA and FL
invsimps_test_df_sorted <- arrange(invsimps_test_df, norep_water_name)
invsimps_test_df_sorted_PA <- filter(invsimps_test_df_sorted, fraction == "Particle")
invsimps_test_df_sorted_FL <- filter(invsimps_test_df_sorted, fraction == "Free")
cor.test(invsimps_test_df_sorted_PA$mean, invsimps_test_df_sorted_FL$mean)

## 
##  Pearson's product-moment correlation
## 
## data:  invsimps_test_df_sorted_PA$mean and invsimps_test_df_sorted_FL$mean
## t = 1.7, df = 10, p-value = 0.12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1378559  0.8235989
## sample estimates:
##       cor 
## 0.4735076

cv_lm_prod_vs_invsimps_PA <- train(
      frac_bacprod ~ mean, data = filter(invsimps, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_prod_vs_invsimps_PA$results       # Particle Community-Wide CV results

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 5.325326 0.7605251 4.051095 2.131925  0.2539022 1.816429

################### Simpson's Evenness ################### 
lm_prod_vs_simpseven_PA <- lm(frac_bacprod ~ mean, data = filter(simpseven, fraction == "Particle"))
cv_lm_prod_vs_simpseven_PA <- train(
      frac_bacprod ~ mean, data = filter(simpseven, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_prod_vs_simpseven_PA$results   # Particle Community-Wide CV results 

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.269156 0.6398939 4.852454 3.012038   0.324493 2.312427

#################################### Per-Capita Production #################################### 
################### Richness ################### 
lm_percell_prod_vs_rich_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, fraction == "Particle"))
summary(lm_percell_prod_vs_rich_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55347 -0.21545 -0.01066  0.12536  0.58830 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.0617648  0.3717671 -21.685 4.44e-09 ***
## mean         0.0023794  0.0006348   3.748  0.00457 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3506 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6095, Adjusted R-squared:  0.5662 
## F-statistic: 14.05 on 1 and 9 DF,  p-value: 0.004567

# Dummy variable regression
summary(lm(log10(fracprod_per_cell_noinf) ~ mean + frac_boolean, 
           data = rich_test_df))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean + frac_boolean, 
##     data = rich_test_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74137 -0.19795 -0.01396  0.24432  0.63068 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -8.3288180  0.2234148 -37.280  < 2e-16 ***
## mean          0.0022273  0.0005832   3.819  0.00107 ** 
## frac_boolean  0.3524489  0.1995645   1.766  0.09264 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3625 on 20 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6973, Adjusted R-squared:  0.667 
## F-statistic: 23.04 on 2 and 20 DF,  p-value: 6.457e-06

PA_rich_percell_residplot <- plot_residuals(lm_percell_prod_vs_rich_PA, filter(richness, fraction == "Particle")$lm(log10(fracprod_per_cell_noinf)), "Log10(Cellular Production) vs Particle-Associated Richness")

## [1] "There are no high leverage points in this model."

## Error in xy.coords(x, y, xlabel, ylabel, log): attempt to apply non-function

cv_lm_percell_prod_vs_rich_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_percell_prod_vs_rich_PA$results      # Particle Per-Capita CV results 

##   intercept      RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
## 1      TRUE 0.4135282 0.6367477 0.3318983 0.1426534   0.327894 0.1262711

################### Shannon Entropy ################### 
lm_percell_prod_vs_shannon_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, fraction == "Particle"))
cv_lm_percell_prod_vs_shannon_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_percell_prod_vs_shannon_PA$results      # Particle Per-Capita CV results 

##   intercept      RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
## 1      TRUE 0.4070696 0.7054441 0.3276488 0.1433208  0.2981486 0.1251932

################### Inverse Simpson ################### 
lm_percell_prod_vs_invsimps_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, fraction == "Particle"))

summary(lm_percell_prod_vs_invsimps_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28651 -0.18384 -0.11125  0.07337  0.56444 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.360961   0.159490 -46.153 5.27e-12 ***
## mean         0.018087   0.003759   4.812 0.000958 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2969 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7201, Adjusted R-squared:  0.689 
## F-statistic: 23.15 on 1 and 9 DF,  p-value: 0.0009581

# Dummy variable regression
summary(lm(log10(fracprod_per_cell_noinf) ~ mean + frac_boolean, data = invsimps_test_df))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean + frac_boolean, 
##     data = invsimps_test_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68100 -0.16952  0.00247  0.13168  0.63767 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -8.01336    0.13664 -58.646  < 2e-16 ***
## mean          0.01819    0.00401   4.536 0.000201 ***
## frac_boolean  0.64869    0.14654   4.427 0.000259 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3347 on 20 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.742,  Adjusted R-squared:  0.7162 
## F-statistic: 28.76 on 2 and 20 DF,  p-value: 1.305e-06

PA_invsimps_percell_residplot <- plot_residuals(lm_percell_prod_vs_invsimps_PA, filter(invsimps, fraction == "Particle")$lm(log10(fracprod_per_cell_noinf)), "Log10(Cellular Production) vs Particle-Associated Inverse Simpson")

## [1] "There are no high leverage points in this model."

## Error in xy.coords(x, y, xlabel, ylabel, log): attempt to apply non-function

cv_lm_percell_prod_vs_invsimps_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_percell_prod_vs_invsimps_PA$results      # Particle Per-Capita CV results 

##   intercept      RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
## 1      TRUE 0.3526888 0.7214097 0.2865713 0.1140385  0.3265452 0.0919178

################### Simpson's Evenness ################### 
lm_percell_prod_vs_simpseven_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, fraction == "Particle"))
cv_lm_percell_prod_vs_simpseven_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))
cv_lm_percell_prod_vs_simpseven_PA$results      # Particle Per-Capita CV results  

##   intercept      RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
## 1      TRUE 0.4368212 0.6842855 0.3526242 0.1383327  0.2979329 0.1103317

Diversity distributions

rich_fraction_wilcox <- wilcox.test(mean ~ fraction, data = richness)

# Make a data frame to draw significance line in boxplot (visually calculated)
rich_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(920, 930, 930, 920)) # WholePart vs WholeFree

rich_vert_distribution_plot <- 
  ggplot(richness, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(150,980), breaks = seq(from = 0, to =925, by = 150)) + 
  ylab("Observed Richness") + xlab("Fraction") + 
  scale_shape_manual(values = season_shapes) +
  geom_path(data = rich_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=980, size = 8, color = "#424645", 
           label= paste("***")) +
  annotate("text", x=1.5, y=880, size = 4, color = "#424645",
           label= paste("p =", round(rich_fraction_wilcox$p.value, digits = 5))) +
  theme(legend.position = "none", axis.title.x = element_blank()) 
  
  
######################################################### INVERSE SIMPSON ######################################################### 
invsimps_fraction_wilcox <- wilcox.test(mean ~ fraction, data = invsimps)

# Make a data frame to draw significance line in boxplot (visually calculated)
invsimps_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(83,85,85,83)) # WholePart vs WholeFree

invsimps_vert_distribution_plot <- 
  ggplot(invsimps, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  geom_jitter(size = 3,  aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  ylab("Inverse Simpson") + xlab("Fraction") + 
  geom_path(data = invsimps_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=80, size = 4, color = "#424645", label= "NS") +
  theme(legend.position = "none", axis.title.x = element_blank()) 

plot_grid(rich_vert_distribution_plot, invsimps_vert_distribution_plot,
          ncol = 2, nrow = 1, labels = c("A", "B")) 

ggsave("div.png", width = 6, height = 3.5)

Unweighted ses.mpd

######################################################### SES MPD DISTRIBUTION: UNWEIGHTED  



################ Figure 3A: Distribution of Particle and Free Unweighted Mean Pairwise distance
## 1. Is there a difference between particle and free Unweighted Mean Pairwise distance?
unweighted_fraction_wilcox <- wilcox.test(mpd.obs.z ~ fraction, data = unweighted_df)
unweighted_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mpd.obs.z by fraction
## W = 30, p-value = 0.01449
## alternative hypothesis: true location shift is not equal to 0

# 2. What are the means of particle and free unweighted PD?
unweighted_df %>%
  group_by(fraction) %>%
  summarize(mean(mpd.obs.z))

## # A tibble: 2 x 2
##   fraction `mean(mpd.obs.z)`
##   <fct>                <dbl>
## 1 Particle            -0.497
## 2 Free                 0.438

## 3. Plot Unweighted Phylogenetic Diversity vs Richness  
# Make a data frame to draw significance line in boxplot (visually calculated)
dat4 <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(1.6,1.7,1.7,1.6)) # WholePart vs WholeFree

unweight_distribution_plot <- 
  ggplot(unweighted_df, aes(y = mpd.obs.z, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_shape_manual(values = season_shapes) + 
  scale_y_continuous(limits = c(-1.5,1.85), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  ylab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  xlab("Fraction") + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 1.5) +
  ##### Particle vs free per-cell production 
  geom_path(data = dat4, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=1.85, size = 8, color = "#424645", label= paste("*"), angle = 90) +
  annotate("text", x=1.5, y=1.35, size = 4, color = "#424645",
           label= paste("p =", round(unweighted_fraction_wilcox$p.value, digits = 2))) +
  theme(legend.position = "none",# axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()


################ Figure 3B: Unweighted Phylogenetic Diversity vs Richness  
## 1. Is there a relationship between richness and Unweighted Mean Pairwise distance?
lm_richness_vs_unweight <- lm(mean ~ mpd.obs.z, data = unweighted_df)

  # Linear model results for particle-associated only 
  summary(lm(mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -191.79 -112.16  -41.89   55.90  278.87 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   515.15      51.16  10.069 1.49e-06 ***
## mpd.obs.z     -83.78      49.91  -1.679    0.124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 155 on 10 degrees of freedom
## Multiple R-squared:  0.2199, Adjusted R-squared:  0.1419 
## F-statistic: 2.818 on 1 and 10 DF,  p-value: 0.1241

  # Linear model results for free-living only 
  summary(lm(mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -100.15  -66.22  -18.70   43.65  172.92 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  337.088     42.753   7.884 1.34e-05 ***
## mpd.obs.z      3.053     77.510   0.039    0.969    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90.18 on 10 degrees of freedom
## Multiple R-squared:  0.0001551,  Adjusted R-squared:  -0.09983 
## F-statistic: 0.001551 on 1 and 10 DF,  p-value: 0.9694

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_richness_vs_unweight<- paste("atop(R^2 ==", round(summary(lm_richness_vs_unweight)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_richness_vs_unweight)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Unweighted Phylogenetic Diversity vs Richness  
unweight_rich_vs_mpd_plot <- 
  ggplot(unweighted_df, aes(y = mean, x = mpd.obs.z)) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3, aes(fill = fraction, shape = season)) + ylab("\n Observed  Richness") +
  scale_shape_manual(values = season_shapes) + 
  xlab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  geom_smooth(method = "lm", color = "#424645", fill = "#424645", alpha = 0.3) + 
  scale_x_continuous(limits = c(-1.5,1.85), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.75, y=750, label=lm_lab_richness_vs_unweight, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position = "none", 
        axis.text.x = element_blank(), axis.title.x = element_blank())


################ Figure 3C: Unweighted Mean Pairwise distance vs Community-Wide Production
# Is there a relationship between Production and Unweighted Mean Pairwise distance?
summary(lm(frac_bacprod ~ mpd.obs.z, data = unweighted_df)) # NS 

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = unweighted_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.416  -9.547  -3.120   8.193  44.030 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.125      3.102   5.520 1.51e-05 ***
## mpd.obs.z      3.901      3.768   1.035    0.312    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.19 on 22 degrees of freedom
## Multiple R-squared:  0.04645,    Adjusted R-squared:  0.003106 
## F-statistic: 1.072 on 1 and 22 DF,  p-value: 0.3118

  # Linear model results for particle-associated only 
  summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.145 -4.312 -1.824  3.402 19.527 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    8.213      2.617   3.138   0.0105 *
## mpd.obs.z     -3.511      2.553  -1.376   0.1990  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.928 on 10 degrees of freedom
## Multiple R-squared:  0.1591, Adjusted R-squared:  0.07501 
## F-statistic: 1.892 on 1 and 10 DF,  p-value: 0.199

  # Linear model results for free-living only 
  summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.996 -14.939   2.192   8.751  37.004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   18.191      8.398   2.166   0.0555 .
## mpd.obs.z     13.410     15.225   0.881   0.3991  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.71 on 10 degrees of freedom
## Multiple R-squared:  0.072,  Adjusted R-squared:  -0.0208 
## F-statistic: 0.7758 on 1 and 10 DF,  p-value: 0.3991

unweight_prod_vs_mpd_plot <- 
  ggplot(unweighted_df, aes(y = frac_bacprod, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3, aes(shape = season)) +
  scale_shape_manual(values = season_shapes) + 
  ylab("\n Community Production \n (μgC/L/day)") +
  xlab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.85), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())



################ Figure 3D: Unweighted Phylogenetic Diversity vs Per Capita (Per-Cell) Production 
## 1. Run the linear model: Is there a relationship between PER-CELL PRODUCTION and Unweighted Mean Pairwise distance?
unweight_vs_percell <- lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = unweighted_df)
summary(unweight_vs_percell)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = unweighted_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72265 -0.28974  0.00388  0.12956  1.31921 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.19878    0.09989 -72.065  < 2e-16 ***
## mpd.obs.z   -0.49728    0.12046  -4.128 0.000478 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4778 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.448,  Adjusted R-squared:  0.4217 
## F-statistic: 17.04 on 1 and 21 DF,  p-value: 0.0004784

summary(aov(log10(fracprod_per_cell_noinf) ~ mpd.obs.z +fraction, data = unweighted_df))

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mpd.obs.z    1  3.890   3.890  20.749 0.000192 ***
## fraction     1  1.044   1.044   5.569 0.028556 *  
## Residuals   20  3.750   0.188                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness

summary(aov(log10(fracprod_per_cell_noinf) ~ fraction + mpd.obs.z, data = unweighted_df))

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## fraction     1  4.139   4.139  22.075 0.000138 ***
## mpd.obs.z    1  0.796   0.796   4.243 0.052673 .  
## Residuals   20  3.750   0.188                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness

lmer(log10(fracprod_per_cell_noinf) ~ mpd.obs.z + (fraction) +
     (0+dummy(fraction, "Particle")), data = Orthodont)

## Error in lmer(log10(fracprod_per_cell_noinf) ~ mpd.obs.z + (fraction) + : could not find function "lmer"

library(lme4)
# Look for the marginal R2 
null<-lmer(log10(fracprod_per_cell_noinf)~(1|fraction),data=unweighted_df)
full<-lmer(log10(fracprod_per_cell_noinf)~mpd.obs.z+(1| fraction),data=unweighted_df)
anova(null,full)

## Data: unweighted_df
## Models:
## null: log10(fracprod_per_cell_noinf) ~ (1 | fraction)
## full: log10(fracprod_per_cell_noinf) ~ mpd.obs.z + (1 | fraction)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## null  3 40.590 43.997 -17.295   34.590                           
## full  4 36.478 41.020 -14.239   28.478 6.1117      1    0.01343 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(full)

## Analysis of Variance Table
##           Df Sum Sq Mean Sq F value
## mpd.obs.z  1  1.092   1.092  5.8239

summary(lm(log10(fracprod_per_cell_noinf) ~ fraction, data = unweighted_df))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ fraction, data = unweighted_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64696 -0.35864  0.04763  0.17090  1.25132 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6.7258     0.1403 -47.946  < 2e-16 ***
## fractionFree  -0.8492     0.1942  -4.373 0.000266 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4652 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4766, Adjusted R-squared:  0.4517 
## F-statistic: 19.12 on 1 and 21 DF,  p-value: 0.0002664

  # Linear model results for particle-associated only 
  summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3795 -0.2419 -0.1651  0.0443  1.1815 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.9248     0.1711 -40.475 1.71e-11 ***
## mpd.obs.z    -0.3299     0.1627  -2.028   0.0732 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4649 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3136, Adjusted R-squared:  0.2373 
## F-statistic: 4.112 on 1 and 9 DF,  p-value: 0.07319

  # Linear model results for free-living only 
  summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63978 -0.29782  0.07694  0.17610  0.76503 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.55622    0.19603 -38.545  3.3e-12 ***
## mpd.obs.z   -0.04303    0.35540  -0.121    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4135 on 10 degrees of freedom
## Multiple R-squared:  0.001463,   Adjusted R-squared:  -0.09839 
## F-statistic: 0.01466 on 1 and 10 DF,  p-value: 0.906

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_percell_unweightPD <- paste("atop(R^2 ==", round(summary(unweight_vs_percell)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(unweight_vs_percell)$coefficients[,4][2]), digits = 4), ")")

## 3. Plot Unweighted Phylogenetic Diversity vs Per Capitra (Per-Cell) Production 
unweight_percell_vs_mpd_plot <- 
  ggplot(unweighted_df, 
       aes(y = log10(fracprod_per_cell_noinf), x = mpd.obs.z)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3,  aes(fill = fraction, shape = season)) + 
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Unweighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) + 
  geom_smooth(method = "lm", color="#424645", fill = "#424645", alpha = 0.3) +
  scale_y_continuous(limits = c(-8.5,-5), breaks = c(-8, -7, -6)) + 
  scale_x_continuous(limits = c(-1.5,1.85), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.75, y=-6, label=lm_lab_percell_unweightPD, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))





######## FIGURE 3
figure3_row1 <- plot_grid(unweight_distribution_plot, unweight_rich_vs_mpd_plot, unweight_prod_vs_mpd_plot, 
          unweight_percell_vs_mpd_plot + theme(legend.position = "none"), 
          labels = c("A", "B", "C", "D"), ncol = 1, nrow = 4,
          rel_heights = c(0.5, 1, 1, 1.2),
          align = "v")

plot_grid(figure3_row1, season_legend,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

Lasso Regressions

unweight_vals <- unweighted_df %>%  
  dplyr::select(mpd.obs.z, norep_filter_name) %>%
  rename(Unweighted_PD = mpd.obs.z)
   
weighted_vals <- weighted_df %>%
  dplyr::select(mpd.obs.z, norep_filter_name) %>%
  rename(Weighted_PD = mpd.obs.z)


wide_all_divs <- otu_alphadiv %>%
  dplyr::select(norep_filter_name ,mean, measure) %>%
  spread(measure, mean)

lasso_data_df <- metadata_pca %>%
  left_join(wide_all_divs, by = "norep_filter_name") %>%
  left_join(unweight_vals, by = "norep_filter_name") %>%
  left_join(weighted_vals, by = "norep_filter_name") %>%
  dplyr::select(-c(project, year, Date, limnion, norep_water_name, dnaconcrep1, perc_attached_bacprod,
                   SD_perc_attached_bacprod, SE_total_bac_abund, SE_perc_attached_abund, SE_attached_bac))

## Warning: Column `norep_filter_name` joining character vector and factor, coercing into character vector

### PARTICLE DATA 
lasso_data_df_particle <- lasso_data_df %>%
  filter(fraction == "Particle" ) %>%
  dplyr::select(-fraction)

lasso_data_df_particle_noprod <- lasso_data_df_particle %>%
  dplyr::select(-c(fracprod_per_cell, fracprod_per_cell_noinf, 
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   norep_filter_name, BGA_cellspermL, Turb_NTU,
                   lakesite, season, station)) 
  
percell_lasso_data_df_particle_noprod <- lasso_data_df_particle %>%
  dplyr::select(-c(fracprod_per_cell, frac_bacprod,
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   lakesite, season, norep_filter_name, station, BGA_cellspermL, Turb_NTU)) %>%
  dplyr::filter(Temp_C != 14.28) # Remove the missing row of data :( )



### FREE DATA 
lasso_data_df_free <- lasso_data_df %>%
  filter(fraction == "Free" ) %>%
  dplyr::select(-fraction)

lasso_data_df_free_noprod <- lasso_data_df_free %>%
  dplyr::select(-c(fracprod_per_cell, fracprod_per_cell_noinf, 
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   norep_filter_name, BGA_cellspermL, Turb_NTU,
                   lakesite, season, station)) 
  
percell_lasso_data_df_free_noprod <- lasso_data_df_free %>%
  dplyr::select(-c(fracprod_per_cell, frac_bacprod,
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   lakesite, season, norep_filter_name, station, BGA_cellspermL, Turb_NTU)) %>%
  dplyr::filter(Temp_C != 14.28) # Remove the missing row of data :( )



## ALL DATA 
all_dat_lasso_comm <- lasso_data_df %>%
    dplyr::select(-c(fracprod_per_cell, fracprod_per_cell_noinf, 
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   norep_filter_name, BGA_cellspermL, Turb_NTU,
                   lakesite, season, station, fraction)) 
  

all_dat_lasso_percapita <- lasso_data_df  %>%
  dplyr::select(-c(fracprod_per_cell, frac_bacprod,
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, fraction,
                   lakesite, season, norep_filter_name, station, BGA_cellspermL, Turb_NTU)) %>%
  dplyr::filter(Temp_C != 14.28) # Remove the missing row of data :( )

# Set the seed for reproducibility of the grid values 
set.seed(777) 

################ PREPARE DATA ################ 
# Subset data needed and scale all o
scaled_comm_data <- 
  lasso_data_df_particle_noprod %>%    # Use data only for the particle samples  
  scale() %>%                          # Scale all of the variables to have mean =0 and sd = 1
  as.data.frame()                      # Make it a dataframe so that model.matrix function works (does not take a matrix)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(frac_bacprod ~ ., scaled_comm_data)[,-1]
y = scaled_comm_data$frac_bacprod
grid = 10^seq(10,-2,length = 100)

################ PREPARE TRAINING & TESTING DATA FOR CROSS VALIDATION ################ 
# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ LASSO ################ 
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min # Minimum lambda
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] TRUE

lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2) # Test

## [1] 1.375897

## Run lasso on the entire dataset with the best lambda value 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -1.478053e-17
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                   .           
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      1.300983e-01
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

## Run lasso on the entire dataset with the best lambda value 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -1.478053e-17
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                   .           
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      1.300983e-01
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

# Set the seed for reproducibility of the grid values 
set.seed(777) 
     
################ PREPARE DATA ################ 
# Subset data needed and scale all o
scaled_comm_data_free <- 
  lasso_data_df_free_noprod %>%        # Use data only for the free-living samples  
  scale() %>%                          # Scale all of the variables to have mean = 0 and sd = 1
  as.data.frame()                      # Make it a dataframe so that model.matrix function works (does not take a matrix)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x <- model.matrix(frac_bacprod ~ ., scaled_comm_data_free)[,-1]
y <- scaled_comm_data_free$frac_bacprod
grid <- 10^seq(10,-2,length = 100)

################ PREPARE TRAINING & TESTING DATA FOR CROSS VALIDATION ################ 
# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ LASSO ################ 
# Run lasso regression with alpha = 1
lasso_divs_train_free_comm <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train_free_comm)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] FALSE

lasso_divs_pred <- predict(lasso_divs_train_free_comm, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2) # Test

## [1] 1.659002

## Run lasso on the entire dataset with the best lambda value 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -1.028226e-15
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                  -5.122124e-01
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           2.183338e-02
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  2.171826e-01
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      .           
## Simpsons_Evenness    .           
## Unweighted_PD        1.397650e-01
## Weighted_PD          .

## Run lasso on the entire dataset with the best lambda value 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -4.292254e-16
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                  -2.218982e-01
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  2.514746e-02
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      .           
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

scaled_percapita_data <- percell_lasso_data_df_particle_noprod %>%   
  dplyr::filter(!is.na(fracprod_per_cell_noinf)) %>%
  mutate(log10_percell = log10(fracprod_per_cell_noinf)) %>%
  dplyr::select(-c(fracprod_per_cell_noinf)) %>%
  scale() %>%
  as.data.frame()
   
set.seed(777)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(log10_percell ~ ., scaled_percapita_data)[,-1]
y = scaled_percapita_data$log10_percell
grid = 10^seq(10,-2,length = 100)

# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ LASSO
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2)) 
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] FALSE

lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 1.743068

## Run lasso on the entire dataset with the best lasso 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -7.322846e-16
## Sample_depth_m       .           
## Temp_C              -1.371144e-01
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                   .           
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund     -5.741093e-03
## attached_bac        -3.957476e-02
## perc_attached_abund  .           
## fraction_bac_abund  -1.164412e-16
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  1.240751e-02
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      4.388114e-01
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

## Run lasso on the entire dataset with the best lasso 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -7.173111e-16
## Sample_depth_m       .           
## Temp_C              -8.717767e-02
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                   .           
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      3.530134e-01
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

scaled_percapita_data_free <- percell_lasso_data_df_free_noprod %>%    
  dplyr::filter(!is.na(fracprod_per_cell_noinf)) %>%
  mutate(log10_percell = log10(fracprod_per_cell_noinf)) %>%
  dplyr::select(-c(fracprod_per_cell_noinf)) %>%
  scale() %>%
  as.data.frame()

set.seed(777)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(log10_percell ~ ., scaled_percapita_data_free)[,-1]
y = scaled_percapita_data_free$log10_percell
grid = 10^seq(10,-2,length = 100)

# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ LASSO
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] TRUE

lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 1.159787

## Run lasso on the entire dataset 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)          1.365298e-15
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                  -1.017130e-01
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      .           
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

## Run lasso on the entire dataset 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se) 

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)          1.365298e-15
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                  -1.017130e-01
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   .           
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      .           
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

# Set the seed for reproducibility of the grid values 
set.seed(777)      
   
################ PREPARE DATA ################ 
# Subset data needed and scale all o
scaled_comm_data_ALL <- 
  all_dat_lasso_comm %>%               # Use community-wide production data only for the all samples  
  scale() %>%                          # Scale all of the variables to have mean =0 and sd = 1
  as.data.frame()                     # Make it a dataframe so that model.matrix function works (does not take a matrix)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x <- model.matrix(frac_bacprod ~ ., scaled_comm_data_ALL)[,-1]
y <- scaled_comm_data_ALL$frac_bacprod
grid <- 10^seq(10,-2,length = 100)

################ PREPARE TRAINING & TESTING DATA FOR CROSS VALIDATION ################ 
# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]


################ LASSO ################ 
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] FALSE

best_lasso_lambda

## [1] 0.3946859

lasso_lambda_1se

## [1] 0.7569721

lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 0.5662862

## Run lasso on the entire dataset with the best lambda value 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                 1
## (Intercept)         -4.328268e-16
## Sample_depth_m       .           
## Temp_C               .           
## SpCond_uSpercm       .           
## TDS_mgperL           .           
## pH                  -2.134918e-01
## ORP_mV               .           
## Chl_Lab_ugperL       .           
## Cl_mgperL            .           
## SO4_mgperL           .           
## NO3_mgperL           .           
## NH3_mgperL           .           
## TKN_mgperL           .           
## SRP_ugperL           .           
## TP_ugperL            .           
## Alk_mgperL           .           
## DO_mgperL            .           
## DO_percent           .           
## total_bac_abund      .           
## attached_bac         .           
## perc_attached_abund  .           
## fraction_bac_abund   4.841677e-02
## PC1                  .           
## PC2                  .           
## PC3                  .           
## PC4                  .           
## PC5                  .           
## PC6                  .           
## Richness             .           
## Shannon_Entropy      .           
## Inverse_Simpson      .           
## Simpsons_Evenness    .           
## Unweighted_PD        .           
## Weighted_PD          .

## Run lasso on the entire dataset with the best lambda value 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se)

## 34 x 1 sparse Matrix of class "dgCMatrix"
##                                1
## (Intercept)         2.081668e-17
## Sample_depth_m      0.000000e+00
## Temp_C              .           
## SpCond_uSpercm      .           
## TDS_mgperL          .           
## pH                  .           
## ORP_mV              .           
## Chl_Lab_ugperL      .           
## Cl_mgperL           .           
## SO4_mgperL          .           
## NO3_mgperL          .           
## NH3_mgperL          .           
## TKN_mgperL          .           
## SRP_ugperL          .           
## TP_ugperL           .           
## Alk_mgperL          .           
## DO_mgperL           .           
## DO_percent          .           
## total_bac_abund     .           
## attached_bac        .           
## perc_attached_abund .           
## fraction_bac_abund  .           
## PC1                 .           
## PC2                 .           
## PC3                 .           
## PC4                 .           
## PC5                 .           
## PC6                 .           
## Richness            .           
## Shannon_Entropy     .           
## Inverse_Simpson     .           
## Simpsons_Evenness   .           
## Unweighted_PD       .           
## Weighted_PD         .

set.seed(777) 
      
scaled_percapita_ALL <- all_dat_lasso_percapita %>%
  dplyr::filter(!is.na(fracprod_per_cell_noinf)) %>%
  mutate(log10_percell = log10(fracprod_per_cell_noinf)) %>%
  dplyr::select(-c(fracprod_per_cell_noinf, perc_attached_abund, fraction_bac_abund, attached_bac)) %>%
  scale() %>%
  as.data.frame()


# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(log10_percell ~ ., scaled_percapita_ALL)[,-1]
y = scaled_percapita_ALL$log10_percell
grid = 10^seq(10,-2,length = 100)

# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train

################ LASSO
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold

plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
# https://stats.stackexchange.com/questions/138569/why-is-lambda-plus-1-standard-error-a-recommended-value-for-lambda-in-an-elastic
lasso_lambda_1se <- cv_lasso_divs$lambda.1se # simplest model whose accuracy is comparable with the best model
best_lasso_lambda == lasso_lambda_1se

## [1] FALSE

lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## Warning in lasso_divs_pred - y_test: longer object length is not a multiple of shorter object length

## Error in mean((lasso_divs_pred - y_test)^2): dims [product 11] do not match the length of object [12]

## Run lasso on the entire dataset 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = best_lasso_lambda, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 31 x 1 sparse Matrix of class "dgCMatrix"
##                               1
## (Intercept)        3.520111e-16
## Sample_depth_m     .           
## Temp_C             .           
## SpCond_uSpercm     .           
## TDS_mgperL         .           
## pH                -6.755466e-02
## ORP_mV             .           
## Chl_Lab_ugperL     .           
## Cl_mgperL          .           
## SO4_mgperL         .           
## NO3_mgperL         .           
## NH3_mgperL         .           
## TKN_mgperL         .           
## SRP_ugperL         .           
## TP_ugperL          .           
## Alk_mgperL         .           
## DO_mgperL          .           
## DO_percent         .           
## total_bac_abund    .           
## PC1                .           
## PC2                .           
## PC3                .           
## PC4                .           
## PC5                .           
## PC6                .           
## Richness           4.154369e-01
## Shannon_Entropy    .           
## Inverse_Simpson    .           
## Simpsons_Evenness  .           
## Unweighted_PD     -6.310163e-02
## Weighted_PD        .

## Run lasso on the entire dataset 
lasso_divs_1se <- glmnet(x, y, alpha = 1, lambda = lasso_lambda_1se, standardize = FALSE)
# What are the lasso coefficients?
predict(lasso_divs_1se, type = "coefficients", s = lasso_lambda_1se)

## 31 x 1 sparse Matrix of class "dgCMatrix"
##                               1
## (Intercept)       -7.854064e-18
## Sample_depth_m     .           
## Temp_C             .           
## SpCond_uSpercm     .           
## TDS_mgperL         .           
## pH                 .           
## ORP_mV             .           
## Chl_Lab_ugperL     .           
## Cl_mgperL          .           
## SO4_mgperL         .           
## NO3_mgperL         .           
## NH3_mgperL         .           
## TKN_mgperL         .           
## SRP_ugperL         .           
## TP_ugperL          .           
## Alk_mgperL         .           
## DO_mgperL          .           
## DO_percent         .           
## total_bac_abund    .           
## PC1                .           
## PC2                .           
## PC3                .           
## PC4                .           
## PC5                .           
## PC6                .           
## Richness           9.842520e-02
## Shannon_Entropy    .           
## Inverse_Simpson    .           
## Simpsons_Evenness  .           
## Unweighted_PD      .           
## Weighted_PD        .

Map of Muskegon Lake

# Load the ggmap package for plotting of Muskegon Lake
library(ggmap)

# Extract the coordinates for Muskegon Lake 
lakesite_coordinates <- read.csv("data/metadata/ML_GPS_Coordinates.csv") %>%
  dplyr::select(lakesite, Latitude, Longitude) %>%
  dplyr::rename(Lat = Latitude, Long = Longitude)

# Set the boundaries for the map 
map <- get_map(c(-86.35, 43.21, -86.235, 43.265),  #left/bottom/right/top
                 zoom = 13, maptype = "toner-background")

# Plot the Muskegon Lake Map 
ggmap(map) + xlab("Longitude") + ylab("Latitude") +
  theme(axis.text = element_text(size = 10, color = "black"),
        axis.title = element_text(size = 12, color = "black", face = "bold")) +
    geom_point(data = lakesite_coordinates, 
               aes(x=Long, y=Lat), color="red", size=4) +
    geom_text(data = lakesite_coordinates,
              aes(x=Long, y=Lat, label=lakesite), col="white", cex=4, vjust = 0.5,hjust = -0.25) 

Faith’s PD

surface_PAFL_otu_pruned_RAREFIED_rm2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1891 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 70 sample variables ]
## tax_table()   Taxonomy Table:    [ 1891 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1891 tips and 1889 internal nodes ]

## Calculate Faith's PD and species richness for sample 
faiths_pd_RAREFIED <- pd(comm_RAREFIED_rm2, phy_RAREFIED_rm2, include.root = FALSE)
faiths_pd_RAREFIED$norep_filter_name <- row.names(faiths_pd_RAREFIED)


# Join Faith's PD with the rest of the metadata 
meta_data_PD <- left_join(faiths_pd_RAREFIED, lasso_data_df, by = "norep_filter_name")

################ Figure S1: species richness and faith's PD
## 1. Run the linear model: Is there a correlation between species richness and faith's PD?
lm_PD_vs_SR <- lm(PD ~ SR, data = meta_data_PD)
summary(lm_PD_vs_SR)

## 
## Call:
## lm(formula = PD ~ SR, data = meta_data_PD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93991 -0.96092  0.03063  1.00013  2.60206 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.836131   0.986599   10.98 2.13e-10 ***
## SR           0.062600   0.002285   27.40  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.52 on 22 degrees of freedom
## Multiple R-squared:  0.9715, Adjusted R-squared:  0.9702 
## F-statistic: 750.9 on 1 and 22 DF,  p-value: < 2.2e-16

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_PD_vs_SR<- paste("atop(R^2 ==", round(summary(lm_PD_vs_SR)$adj.r.squared, digits = 3), ",",
             "p ==", round(unname(summary(lm_PD_vs_SR)$coefficients[,4][2]), digits = 20), ")")

## 3. Plot Figure S1: species richness and faith's PD
ggplot(meta_data_PD, aes(y = PD, x = SR)) + 
  geom_point(size = 3) + 
  ylab("Faith's PD") + 
  xlab("Species Richness") +
  geom_smooth(method = "lm", color = "black") +
  # Add the R2 and p-value to the plot 
  annotate("text", x=600, y=35, label=lm_lab_PD_vs_SR, parse = TRUE, color = "black", size = 4) 

Correlations between Particle and Free of values in Figure 1

work_df <- metadata %>%
  dplyr::select(norep_filter_name, fraction, fraction_bac_abund, frac_bacprod, fracprod_per_cell_noinf) %>%
  mutate(norep_water_name = paste(substr(norep_filter_name, 1,4), substr(norep_filter_name, 6,9), sep = "")) %>%
  dplyr::select(-norep_filter_name)

part_work_df <- work_df %>%
  filter(fraction == "Particle") %>%
  rename(part_bacabund = fraction_bac_abund,
         part_prod = frac_bacprod, 
         part_percell_prod = fracprod_per_cell_noinf) %>%
  dplyr::select(-fraction)

free_work_df <- work_df %>%
  filter(fraction == "Free") %>%
  rename(free_bacabund = fraction_bac_abund,
         free_prod = frac_bacprod, 
         free_percell_prod = fracprod_per_cell_noinf) %>%
  dplyr::select(-fraction)

byfrac_df <- part_work_df %>%
  left_join(free_work_df, by = "norep_water_name")

byfrac_df$season <- substr(byfrac_df$norep_water_name, 5,5) # 7th letter = month sampled
byfrac_df$season <- as.character(byfrac_df$season)
byfrac_df$season <- ifelse(byfrac_df$season == "5", "Spring", 
                             ifelse(byfrac_df$season == "7", "Summer", 
                                    ifelse(byfrac_df$season == "9", "Fall",
                                           "NA")))
byfrac_df$season <- factor(byfrac_df$season, levels = c("Spring", "Summer", "Fall"))

summary(lm(log10(part_bacabund) ~ log10(free_bacabund), data = filter(byfrac_df, norep_water_name != "MOTE515")))

## 
## Call:
## lm(formula = log10(part_bacabund) ~ log10(free_bacabund), data = filter(byfrac_df, 
##     norep_water_name != "MOTE515"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52690 -0.08048  0.05903  0.19565  0.35255 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)            2.0900     3.2247   0.648    0.533
## log10(free_bacabund)   0.4264     0.5532   0.771    0.461
## 
## Residual standard error: 0.3111 on 9 degrees of freedom
## Multiple R-squared:  0.06194,    Adjusted R-squared:  -0.04229 
## F-statistic: 0.5943 on 1 and 9 DF,  p-value: 0.4605

plot1 <- ggplot(filter(byfrac_df, norep_water_name != "MOTE515"), 
       aes(x = log10(free_bacabund), y = log10(part_bacabund))) +
  xlab("Free") +  ylab("Particle") + 
  ggtitle("Log10(Bacterial Counts) \n (cells/mL)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  theme(legend.position = "bottom",
        legend.title = element_blank())




################ Community-Wide production correlation between particle and free
## 1. Run the linear regression 
lm_prod_corr <- lm(part_prod ~ free_prod, data = byfrac_df)

# 2. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_PAvsFL <- paste("atop(R^2 ==", round(summary(lm_prod_corr)$adj.r.squared, digits = 3), ",",
            "p ==", round(unname(summary(lm_prod_corr)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Community production correlation between particle and free
plot2 <- ggplot(byfrac_df, aes(x = free_prod, y = part_prod)) +
  xlab("Free") +  ylab("Particle") + 
  ggtitle("Community Production \n(μgC/L/day)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  geom_smooth(method = "lm", color = "black")  + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=20, y=30, label=lm_lab_perliter_PAvsFL, parse = TRUE, color = "black", size = 4) +
  theme(legend.position = "none")


################ Per-Capita production correlation between particle and free
## 1. Run the linear regression 
lm_percell_corr <- lm(log10(part_percell_prod) ~ log10(free_percell_prod), data = byfrac_df)

# 2. Extract the R2 and p-value from the linear model: 
lm_lab_percell_PAvsFL <- paste("atop(R^2 ==", round(summary(lm_percell_corr)$adj.r.squared, digits = 3), ",",
            "p ==", round(unname(summary(lm_percell_corr)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Per-Capita production correlation between particle and free
plot3 <- ggplot(byfrac_df,
       aes(x = log10(free_percell_prod), y = log10(part_percell_prod))) +  
  xlab("Free") +  ylab("Particle") + 
  ggtitle("log10(Per-Capita Production)\n (μgC/cell/day)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  geom_smooth(method = "lm", color = "black") + 
  # Add the R2 and p-value to the plot 
  annotate("text", y = -5.8, x=-7.9, label=lm_lab_percell_PAvsFL, parse = TRUE, color = "black", size = 4) +
  theme(legend.position = "none")


# Make the final multi-paneled Plot
# Extract the legend
legend_s1 <- get_legend(plot1)

# Make the 3 plots
top_row_S1 <- plot_grid(plot1 +theme(legend.position = "none"), plot2, plot3, 
          nrow = 1, ncol = 3, 
          labels = c("A", "B", "C"), 
          align = "h")

# Put the legend and 3 plots together
plot_grid(top_row_S1, legend_s1,
          rel_heights = c(1, 0.05),
          nrow = 2, ncol = 1)

Shannon and Simpson’s Evenness

######################################################### SHANNON ENTROPY ######################################################### 
shannon_fraction_wilcox <- wilcox.test(mean ~ fraction, data = shannon)
shannon_fraction_wilcox   

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 119, p-value = 0.00556
## alternative hypothesis: true location shift is not equal to 0

filter(shannon) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##   <fct>           <dbl>      <dbl>
## 1 Particle         4.53      0.559
## 2 Free             3.98      0.296

# Make a data frame to draw significance line in boxplot (visually calculated)
shannon_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(5.8, 5.9, 5.9, 5.8)) # WholePart vs WholeFree

shannon_distribution_plot <- 
  ggplot(shannon, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(shape = season, fill = fraction), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  xlab("Shannon Entropy") +
  xlab("Fraction") + 
  scale_shape_manual(values = season_shapes) +
    # Add line and pval
  geom_path(data = shannon_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=6.1, size = 8, color = "#424645", label= paste("**"), angle = 90) +
  annotate("text", x=1.5, y=5.5, size = 4, color = "#424645",
           label= paste("p =", round(shannon_fraction_wilcox$p.value, digits = 3))) +
  theme(legend.position = "none",# axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()

################ Shannon vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_shannon_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_shannon_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_shannon_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Shannon vs Community-wide (Per-Liter) Production 
prod_vs_shannon_plot <-  
  ggplot(shannon, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n Community Production \n (μgC/L/day)") + 
  scale_x_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  xlab("Shannon Entropy") +
  geom_smooth(data=filter(shannon, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=45, label=lm_lab_perliter_shannon_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 


################ Shannon vs Per Capitra (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_shannon_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_shannon_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_shannon_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Shannon vs Per Capitra (Per-Cell) Production 
percell_prod_vs_shannon_plot <- 
  ggplot(shannon, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Shannon Entropy") +
  geom_smooth(data=filter(shannon, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  #scale_y_continuous(limits = c(-8e-7,5e-6), breaks = seq(from = 0, to = 6e-7, by = 3e-7)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=-7.5, label=lm_lab_percell_shannon_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position = "none")



shannon_plots <- plot_grid(shannon_distribution_plot, prod_vs_shannon_plot, percell_prod_vs_shannon_plot,
          labels = c("A", "C", "E"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.85, 1.1),
          align = "v")


######################################################### SIMPSON'S EVENNESS ######################################################### 
simpseven_fraction_wilcox <- wilcox.test(mean ~ fraction, data = simpseven)
simpseven_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 55, p-value = 0.3474
## alternative hypothesis: true location shift is not equal to 0

filter(simpseven) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##   <fct>           <dbl>      <dbl>
## 1 Particle       0.0589     0.0254
## 2 Free           0.0714     0.0172

# Make a data frame to draw significance line in boxplot (visually calculated)
simpseven_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(0.14, 0.15, 0.15, 0.14)) # WholePart vs WholeFree

simpseven_distribution_plot <- 
  ggplot(simpseven, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  ylab("Simpson's Evenness") +
  scale_shape_manual(values = season_shapes) +
  xlab("Fraction") + 
  geom_path(data = simpseven_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=0.14, size = 4, color = "#424645", label= "NS") +
  theme(legend.position = "none",
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()




################ Simpson's Evenness vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_simpseven_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_simpseven_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_simpseven_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Simpson's Evenness vs Community-wide (Per-Liter) Production 
prod_vs_simpseven_plot <-  
  ggplot(simpseven, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n Community Production \n (μgC/L/day)") + 
  ylab("Simpson's Evenness") +
  geom_smooth(data=filter(simpseven, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.13, y=15, label=lm_lab_perliter_simpseven_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 



################ Simpson's Evenness vs Per Capitra (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_simpseven_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_simpseven_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_simpseven_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Simpson's Evenness vs Per Capitra (Per-Cell) Production 
percell_prod_vs_simpseven_plot <- 
  ggplot(simpseven, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Simpson's Evenness") +
  geom_smooth(data=filter(simpseven, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.13, y=-6.3, label=lm_lab_percell_simpseven_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position ="none")


simpseven_plots <- plot_grid(simpseven_distribution_plot, prod_vs_simpseven_plot, percell_prod_vs_simpseven_plot,
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.85, 1.1),
          align = "v")


#### Remove the axis labels for prettier multipaneled plotting
simpseven_plots_noyaxis <- 
  plot_grid(simpseven_distribution_plot + theme(axis.title.y = element_blank(), axis.text.y = element_blank()), 
            prod_vs_simpseven_plot + theme(axis.title.y = element_blank()), 
            percell_prod_vs_simpseven_plot + theme(axis.title.y = element_blank()), 
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.8, 1.1),
          align = "v")

figS2_row1 <- plot_grid(shannon_plots, simpseven_plots_noyaxis,
          ncol = 2, nrow = 1, rel_widths = c(1, 0.75),
          align = "h")

###### PLOT FIGURE S5
plot_grid(figS2_row1, season_legend,
          ncol = 1, nrow = 2, rel_heights = c(1, 0.05))

Compare Diversity Estimates across Station and Season

long_div_df <- lasso_data_df %>%  
  dplyr::select(norep_filter_name, lakesite:season, Richness:Weighted_PD) %>%
  gather(measure, mean, Richness:Weighted_PD) %>%
  mutate(measure = factor(measure, levels = 
                            c("Richness", "Shannon_Entropy", "Inverse_Simpson", "Simpsons_Evenness", "Unweighted_PD", "Weighted_PD")))


divs_PAFLA_plot_lakesite <- ggplot(long_div_df, aes(y = mean, x = lakesite, shape = lakesite)) +
  ylab("Mean Diversity Value\n By Lake Station") +
  facet_wrap(~measure, scale = "free_y", ncol = 6) +
    geom_point(size = 3, aes(fill = fraction), color = "black", position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "black", aes(fill = fraction)) +
  scale_fill_manual(values = fraction_colors) + 
  scale_shape_manual(values = lakesite_shapes) + 
  scale_color_manual(values = fraction_colors) + 
  theme(legend.position = "none", axis.title.x = element_blank(), 
        axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1)) 

plot_A <- plot_grid(divs_PAFLA_plot_lakesite, lakesite_legend,
          nrow = 2, ncol = 1,  labels = c("A",""),
          rel_heights = c(1, 0.05))


divs_PAFLA_plot_season <- ggplot(long_div_df, aes(y = mean, x = season, shape = season)) +
  ylab("Mean Diversity Value \n By Season") +
  facet_wrap(~measure, scale = "free_y", ncol = 6) +
   geom_point(size = 3, aes(fill = fraction), color = "black", position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "black", aes(fill = fraction)) +
  scale_fill_manual(values = fraction_colors) + 
  scale_shape_manual(values = season_shapes) + 
  scale_color_manual(values = fraction_colors) + 
  theme(legend.position = "none", axis.title.x = element_blank(), 
        axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1)) 

plot_B <- plot_grid(divs_PAFLA_plot_season, season_legend,
          nrow = 2, ncol = 1, labels = c("B",""),
          rel_heights = c(1, 0.05))

###### PLOT FIGURE S6
plot_grid(plot_A, plot_B,
          nrow = 2, ncol = 1)

See analysis/OTU_Removal_Analysis* files on the GitHub Repo

Per capita production combined models

################ Figure S7A: Per-cell Heterotrophic Production vs All Observed Richness
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_rich_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = richness))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = richness))$coefficients[,4][2]), digits = 6), ")")

## 2.  Figure S7A: Per-cell Heterotrophic Production vs All Observed Richness
plot_all_rich_percell <- 
  ggplot(richness, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") + xlab("Observed Richness") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(150,950), breaks = seq(from = 0, to =925, by = 150)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=790, y=-8, label=lm_lab_overall_rich_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))


################ Figure S7B: Per-cell Heterotrophic Production vs All shannon
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_shannon_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = shannon))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = shannon))$coefficients[,4][2]), digits = 5), ")")

## 2.  Figure S7B: Per-cell Heterotrophic Production vs All Shannon
plot_all_shannon_percell <- 
  ggplot(shannon, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Shannon Entropy") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(3.5,6), breaks = seq(from = 3, to =6, by = 0.5)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=-8, label=lm_lab_overall_shannon_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))

################ Figure S7C: Per-cell Heterotrophic Production vs All Inverse Simpson
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_invsimps_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = invsimps))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = invsimps))$coefficients[,4][2]), digits = 5), ")")

## 2.  Figure S7C: Per-cell Heterotrophic Production vs All Inverse Simpson
plot_all_invsimps_percell <- 
  ggplot(invsimps, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Inverse Simpson") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=70, y=-8, label=lm_lab_overall_invsimps_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))

################ Figure S7D: Per-cell Heterotrophic Production vs All Simpson's Evenness
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_simpseven_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = simpseven))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = simpseven))$coefficients[,4][2]), digits = 2), ")")

## 2.  Figure S7D: Per-cell Heterotrophic Production vs All Inverse Simpson
plot_all_simpseven_percell <- 
  ggplot(simpseven, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Simpson's Evenness") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.12, y=-8, label=lm_lab_overall_simpseven_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))


figs3_row1 <- plot_grid(plot_all_rich_percell + theme(legend.position = "none"), 
          plot_all_shannon_percell + theme(axis.title.y = element_blank(), legend.position = "none"), 
          plot_all_invsimps_percell + theme(axis.title.y = element_blank(), legend.position = "none"), 
          plot_all_simpseven_percell + theme(axis.title.y = element_blank(), legend.position = "none"),
          align = "h", labels = c("A", "B", "C", "D"),
          rel_widths = c(1.05, 0.9, 0.9, 0.9),
          nrow = 1, ncol = 4)

######## FIGURE S7
# legend
legend_grey <- get_legend(plot_all_rich_percell)


plot_grid(figs3_row1, legend_grey,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

Weighted ses.mpd

######################################################### SES MPD DISTRIBUTION: WEIGHTED 
weighted_fraction_wilcox <- wilcox.test(mpd.obs.z ~ fraction, data = weighted_df)
weighted_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mpd.obs.z by fraction
## W = 81, p-value = 0.6297
## alternative hypothesis: true location shift is not equal to 0

filter(weighted_df) %>%
  group_by(fraction) %>%
  summarize(mean(mpd.obs.z))

## # A tibble: 2 x 2
##   fraction `mean(mpd.obs.z)`
##   <fct>                <dbl>
## 1 Particle            -0.361
## 2 Free                -0.385

# Make a data frame to draw significance line in boxplot (visually calculated)
dat5 <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(1.6,1.7,1.7,1.6)) # WholePart vs WholeFree

weight_distribution_plot <- 
  ggplot(weighted_df, aes(y = mpd.obs.z, x = fraction)) +
  #scale_color_manual(values = fraction_colors) + 
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3.5, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_shape_manual(values = season_shapes) + 
  scale_y_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  ylab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  xlab("Fraction") + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 1.5) +
  ##### Particle vs free per-cell production 
  geom_path(data = dat5, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=1.55, fontface = "bold",  size = 3.5, color = "#424645", label= "NS") +
  theme(legend.position = "none", #axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()


# Is there a relationship between inverse simpson and Weighted Mean Pairwise distance?
weight_invsimps_vs_mpd_plot <- 
  ggplot(weighted_df, aes(y = mean, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + ylab("Inverse Simpson") +
  scale_shape_manual(values = season_shapes) + 
  #xlab("Standardized Effect Size \n Weighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())




# Is there a relationship between PER-LITER PRODUCTION and WEIGHTED Mean Pairwise distance?
weight_prod_vs_mpd_plot <- 
  ggplot(weighted_df, aes(y = frac_bacprod, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + 
  scale_shape_manual(values = season_shapes) + 
  ylab("\n Community Production \n (μgC/L/day)") + 
  xlab("Standardized Effect Size \n Weighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  scale_color_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())



################ Weighted Phylogenetic Diversity vs Per Capita (Per-Cell) Production 
weight_percell_vs_mpd_plot <- 
  ggplot(weighted_df, 
       aes(y = log10(fracprod_per_cell_noinf), x = mpd.obs.z, fill = fraction)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + 
  scale_shape_manual(values = season_shapes) + 
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Weighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  scale_y_continuous(limits = c(-8.5,-5), breaks = c(-8, -7, -6)) + 
  theme(legend.title = element_blank(), legend.position ="none", 
        legend.text = element_text(size = 14))


row1_weighted_PD_plot <- plot_grid(weight_distribution_plot, weight_invsimps_vs_mpd_plot, weight_prod_vs_mpd_plot, weight_percell_vs_mpd_plot, 
          labels = c("A", "B", "C", "D"), ncol = 1, nrow = 4,
          rel_heights = c(0.5, 0.8, 0.8, 1.2),
          align = "v")

### PLOT FIGURE S8
plot_grid(row1_weighted_PD_plot, season_legend,
          ncol = 1, nrow = 2,
          rel_heights = c(1, 0.05))

Randomized Richness

#########################################################  Subset only richness data 
### These are the richness values for the fake samples 
#rich_stats <- filter(otu_alphadiv, measure == "Richness") %>%
#  dplyr::select(1:2) %>%
#  rename(mean_richness = mean) %>%
#  mutate(sample = paste("Sample_", seq(1:nrow(filter(otu_alphadiv, measure == "Richness"))), sep = ""),
#         mean_richness = matround(mean_richness))

## Pick OTUs to match these richness values 

  # List the otus from ALL samples 
#  all_otus <- taxa_names(surface_PAFL_otu_pruned_rm2)
  
  # Obtain the OTU table from the phyloseq object
#  otutab <- otu_table(surface_PAFL_otu_pruned_rm2)
  # Make all the counts to be 0
#  otutab_newvals <- apply(otutab, c(1, 2), function(x) 0)

  # Stop if things are wrong 
#  stopifnot(colnames(otutab_newvals) == all_otus)                       # Do the OTU names match?
#  stopifnot(rownames(otutab_newvals) == rich_stats$norep_filter_name)   # Do the sample names match?
  
# Make it reproducible!   
#set.seed(777)
  
#for (row in 1:nrow(rich_stats)) {
  
  # Pick the richness value 
#  rich_val <- rich_stats[row, 2]  
  
  # Sample the OTUs to represent that richness value 
#  col_index <- sample(ncol(otutab_newvals), rich_val, replace = FALSE, prob = NULL)
  
  # make all other columns 0
#  otutab_newvals[row, col_index] <- 1

#}


## Calculate the tree for those randomized samples 
# create a new phyloseq object 
#random_physeq_presab_raw <- phyloseq(otu_table(otutab_newvals, taxa_are_rows = FALSE), 
#                                 tax_table(surface_PAFL_otu_pruned_rm2), sample_data(surface_PAFL_otu_pruned_rm2))
#random_physeq_presab_raw

# Remove taxa that are 0!
#random_physeq_presab_pruned <- prune_taxa(taxa_sums(random_physeq_presab_raw) > 0, random_physeq_presab_raw) 
#random_physeq_presab_pruned

# Calculate tree 
# Obtain the OTU names that were retained in the dataset
#tax <- data.frame(tax_table(random_physeq_presab_pruned))
#vector_of_otus <- as.vector(tax$OTU)

# Write out the file for processing/fasttree
# write(vector_of_otus, file = "data/PhyloTree/randomized/random_physeq_presab_pruned_OTUnames.txt", ncolumns = 1, append = FALSE, sep = "\n")

# Read in the tree
randomized_tree <- read.tree(file = "data/PhyloTree/randomized/newick_tree_randomized_rmN.tre")
  
  #random_physeq_presab_pruned_tree <- merge_phyloseq(random_physeq_presab_pruned, randomized_tree)
  #random_physeq_presab_pruned_tree
#save(list="random_physeq_presab_pruned_tree", file=paste0("data/PhyloTree/randomized/random_physeq_presab_pruned_tree.RData")) 

load("data/PhyloTree/randomized/random_physeq_presab_pruned_tree.RData")
random_physeq_presab_pruned_tree

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2911 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 44 sample variables ]
## tax_table()   Taxonomy Table:    [ 2911 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2911 tips and 2909 internal nodes ]

# First force the OTU 
randomized_otu <- matrix(otu_table(random_physeq_presab_pruned_tree), nrow = nrow(otu_table(random_physeq_presab_pruned_tree)))
rownames(randomized_otu) <- sample_names(random_physeq_presab_pruned_tree)
colnames(randomized_otu) <- taxa_names(random_physeq_presab_pruned_tree)
    
  
## Calculate input for SES_MPD  
# Convert the presence/absence data to standardized abundanced  vegan function `decostand' , NOTE: method = "pa"
otu_decostand <- decostand(randomized_otu, method = "pa")
# check total abundance in each sample
apply(otu_decostand, 1, sum)

## MBREJ515 MBREJ715 MBREJ915 MBREK515 MBREK715 MBREK915 MDPEJ515 MDPEJ715 MDPEJ915 MDPEK515 MDPEK715 MDPEK915 MINEJ515 MINEJ715 MINEJ915 MINEK515 MINEK715 MINEK915 MOTEJ515 MOTEJ715 MOTEJ915 MOTEK515 MOTEK715 MOTEK915 
##      906      574      434      268      256      358      493      447      476      276      284      381      840      632      586      452      383      511      505      343      444      274      238      381

# check for mismatches/missing species between community data and phylo tree
randomized_matches <- match.phylo.comm(randomized_tree, otu_decostand)
# the resulting object is a list with $phy and $comm elements.  replace our
# original data with the sorted/matched data
phy_randomized_rm2 <- randomized_matches$phy
comm_randomized_rm2 <- randomized_matches$comm

# Calculate the phylogenetic distances
phy_dist_randomized_rm2 <- cophenetic(phy_randomized_rm2)

  
## Calculate SES_MPD
###################################### INDEPENDENT SWAP ############################################
# calculate standardized effect size mean pairwise distance (ses.mpd)
unweighted_sesMPD_indepswap_randomized <- ses.mpd(comm_randomized_rm2, phy_dist_randomized_rm2, null.model = "independentswap", 
                                     abundance.weighted = FALSE, runs = 999)

df <- unweighted_sesMPD_indepswap_randomized

df$names <- row.names(df)
df_2 <- make_metadata_norep(df) %>%
  mutate(fraction = fct_recode(fraction, Particle = "WholePart", Free = "WholeFree"))
  
  
# Is there a relationship between richness and Unweighted Mean Pairwise distance?
summary(lm(ntaxa ~ mpd.obs.z, data = df_2))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = df_2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -214.61 -128.21  -13.53   56.63  465.39 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   448.65      35.29  12.714  1.3e-11 ***
## mpd.obs.z      75.24      95.94   0.784    0.441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 172.8 on 22 degrees of freedom
## Multiple R-squared:  0.02719,    Adjusted R-squared:  -0.01703 
## F-statistic: 0.615 on 1 and 22 DF,  p-value: 0.4413

summary(lm(ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == "Particle")))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -208.89 -109.98  -51.44   39.09  346.60 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   557.41      50.97  10.935 6.96e-07 ***
## mpd.obs.z     -18.67     153.80  -0.121    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 175.3 on 10 degrees of freedom
## Multiple R-squared:  0.001471,   Adjusted R-squared:  -0.09838 
## F-statistic: 0.01474 on 1 and 10 DF,  p-value: 0.9058

summary(lm(ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == "Free")))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -108.70  -51.56  -17.29   57.16  160.49 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   343.12      25.00  13.727 8.18e-08 ***
## mpd.obs.z      68.10      62.35   1.092      0.3    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 85.34 on 10 degrees of freedom
## Multiple R-squared:  0.1066, Adjusted R-squared:  0.01726 
## F-statistic: 1.193 on 1 and 10 DF,  p-value: 0.3003

randomized_rich_plot <- ggplot(df_2, aes(y = ntaxa, x = mpd.obs.z, fill = fraction)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + ylab("Randomized Richness") +
  xlab("Unweighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) + 
  scale_x_continuous(limits = c(-1,1)) + 
  theme(legend.title = element_blank(), legend.position = "none")


plot_grid(randomized_rich_plot, season_legend,
          rel_heights = c(1, 0.05), nrow = 2, ncol = 1)

PCA with Euclidian Distances

# Make the Supplemental Figure
par(mar = c(4.5, 4.5, 1, 1))
biplot(pca_environ,
     xlab = paste("PC1", paste(round(summary(pca_environ)$cont$importance[2,1]*100, digits = 2), "%", sep = ""), sep = ": "),
     ylab = paste("PC2", paste(round(summary(pca_environ)$cont$importance[2,2]*100, digits = 2), "%", sep = ""), sep = ": "))

par(mar = c(5,5,2,5))
plot(summary(pca_environ)$cont$importance[2,]*100, 
     xlab = "PCA Axis", 
     ylab = "Variation Explained Per Axis",
     ylim = c(0, 105),
     col = "cornflowerblue",
     cex =2,
     pch = 16)
par(new = T)
plot(summary(pca_environ)$cont$importance[3,]*100,
       cex = 2,
       pch = 17,
       ylim = c(0, 105),
       col = "firebrick3",
       axes=F, 
       xlab=NA, 
       ylab=NA)
axis(side = 4)
mtext(side = 4, line = 3, "Total Accumulated Variation")
legend("right",
       legend=c("Per Axis", "Total"),
       pch=c(16, 17), col=c("cornflowerblue", "firebrick3"))

# Amount of variation explained by the first two axes
summary(pca_environ)$cont$importance[3, 2]*100

## [1] 69.271

The first two axes of the PCA explain ~70% of the variation in the environmental data!

Prepare data frames for table S1 and S2

# Add the phylogenetic diversity to the dataframes
weighted_PD_df <- weighted_df %>%
  dplyr::select(norep_filter_name, mpd.obs.z) %>%
  rename(weighted_PD = mpd.obs.z)

unweighted_PD_df <- unweighted_df %>%
  dplyr::select(norep_filter_name, mpd.obs.z) %>%
  rename(unweighted_PD = mpd.obs.z)

metadata_pca_PD <- metadata_pca %>%
  left_join(unweighted_PD_df, by = "norep_filter_name") %>%
  left_join(weighted_PD_df, by = "norep_filter_name")
  

## Prepare the dataframes for Free and Particle data
metadata_pca_PD_free <- filter(metadata_pca_PD, fraction == "Free")
metadata_pca_PD_part <- filter(metadata_pca_PD, fraction == "Particle")  

# List of environmental variables to run linear models on 
        # BGA_cellspermL and Turb_NTU do not have enough data points for regression
lm_variables <- c("frac_bacprod", "fracprod_per_cell_noinf", "Temp_C", "SpCond_uSpercm", "TDS_mgperL", 
                  "pH", "ORP_mV", "Chl_Lab_ugperL", "Cl_mgperL", "SO4_mgperL", "NO3_mgperL", 
                  "NH3_mgperL", "TKN_mgperL", "SRP_ugperL", "TP_ugperL", "Alk_mgperL", "DO_mgperL", 
                  "total_bac_abund","attached_bac", "dnaconcrep1", 
                  "PC1", "PC2", "PC3" ,"PC4", "PC5", "PC6",
                  "unweighted_PD", "weighted_PD") 

########################################################
##############  All Samples: Community-wide 
lm_summary_all_comm <- metadata_pca_PD %>%
  dplyr::select(one_of(lm_variables), -fracprod_per_cell_noinf) %>%
  gather(key = independent, value = measurement, -frac_bacprod) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(frac_bacprod ~ measurement, data = .))) %>% 
  mutate(fraction = "All Samples", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  All Samples: Per-capita  
lm_summary_all_percap <- metadata_pca_PD %>%
  dplyr::select(one_of(lm_variables), -frac_bacprod) %>%
  filter(!is.na(fracprod_per_cell_noinf)) %>%   # Remove samples that are not NA
  gather(key = independent, value = measurement, -fracprod_per_cell_noinf) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ measurement, data = .))) %>% 
  mutate(fraction = "All Samples", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 



########################################################
##############  Particle: Community-wide 
lm_summary_particle_comm <- metadata_pca_PD_part %>%
  dplyr::select(one_of(lm_variables), -fracprod_per_cell_noinf) %>%
  gather(key = independent, value = measurement, -frac_bacprod) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(frac_bacprod ~ measurement, data = .))) %>% 
  mutate(fraction = "Particle", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  Particle: Per-capita  
lm_summary_particle_percap <- metadata_pca_PD_part %>%
  dplyr::select(one_of(lm_variables), -frac_bacprod) %>%
  filter(!is.na(fracprod_per_cell_noinf)) %>%   # Remove samples that are not NA
  gather(key = independent, value = measurement, -fracprod_per_cell_noinf) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ measurement, data = .))) %>% 
  mutate(fraction = "Particle", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 


########################################################
##############  Free: Community-wide 
lm_summary_free_comm <- metadata_pca_PD_free %>%
  dplyr::select(one_of(lm_variables), -fracprod_per_cell_noinf) %>%
  gather(key = independent, value = measurement, -frac_bacprod) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(frac_bacprod ~ measurement, data = .))) %>% 
  mutate(fraction = "Free", dependent = "Community Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 

##############  Free: Per-capita  
lm_summary_free_percap <- metadata_pca_PD_free %>%
  dplyr::select(one_of(lm_variables), -frac_bacprod) %>%
  filter(!is.na(fracprod_per_cell_noinf)) %>%   # Remove samples that are not NA
  gather(key = independent, value = measurement, -fracprod_per_cell_noinf) %>%
  mutate(measurement = as.numeric(measurement)) %>%
  group_by(independent) %>% 
  do(glance(lm(log10(fracprod_per_cell_noinf) ~ measurement, data = .))) %>% 
  mutate(fraction = "Free", dependent = "Per-Capita Production") %>%
  dplyr::select(fraction, dependent, independent, AIC, adj.r.squared, p.value) %>%
  ungroup() %>%
  mutate(FDR.p = p.adjust(p.value, method = "fdr")) 


## Filtered PCA Linear Model Results
# Put all of the PCA and environmental  linear model results together into one dataframe 
all_pca_environ_lm_results <- 
  bind_rows(lm_summary_particle_comm, lm_summary_particle_percap,
          lm_summary_free_comm, lm_summary_free_percap,
          lm_summary_all_comm, lm_summary_all_percap) %>%
  dplyr::rename(independent_var = independent,
                dependent_var = dependent) 

## Combine the results with the diversity linear models 
individual_lm_results <- bind_rows(lm_div_results, all_pca_environ_lm_results) %>%
  mutate(AIC = round(AIC, digits = 2),  
         adj.r.squared = round(adj.r.squared, digits = 2),
         p.value = round(p.value, digits = 4),
         FDR.p = round(FDR.p, digits = 4)) %>%
  rename(FDR.p.value = FDR.p)

All Regression Results

### Per-Capita Production     
datatable(individual_lm_results, 
       options = list(pageLength = 20), rownames = FALSE)